rm(list = ls())
Sys.info()[7]

N_bootstrap <- 3 # batch size
BATCH <- 1

# Set file paths 
DATA <- "../external_links/real"
TEMP <- "../temp"
RESULTS <- "../output"
MODEL_RESULTS <- "../external/estimate_w_xs" # for separate dir

# Go to code folder
setwd("./")
sink(paste0("bootstrap_se_compile.txt"))

folder <- paste0(TEMP, "/bs_samples")
if (file.exists(folder)) {
  cat("The folder already exists")
} else {
  dir.create(folder)
}

Sys.info()[7]
Sys.time()

# Packages
# library(MASS)
# library(tmvtnorm)
# library(haven)
# library(dplyr)
# library(tidyr)
# library(nloptr)
# library(corpcor)
# library(stargazer)
# library(pracma)
# library(foreign)
# library(writexl)
# library(broom)
# library(png)
# library(tidyverse)

print("Begin time")
start_time <- Sys.time()
print(start_time)

real_bs_table_1  <- read.csv(paste0(RESULTS, "/bootstrap_se_batch1.csv"), header = FALSE)
real_bs_table_2  <- read.csv(paste0(RESULTS, "/bootstrap_se_batch2.csv"), header = FALSE)
real_bs_table_3  <- read.csv(paste0(RESULTS, "/bootstrap_se_batch3.csv"), header = FALSE)
real_bs_table_4  <- read.csv(paste0(RESULTS, "/bootstrap_se_batch4.csv"), header = FALSE)
real_bs_table_5  <- read.csv(paste0(RESULTS, "/bootstrap_se_batch5.csv"), header = FALSE)

real_bs_table <- rbind(real_bs_table_1, real_bs_table_2, real_bs_table_3, real_bs_table_4, real_bs_table_5)
# real_bs_table <- real_bs_table[,-1]
sd_mu_a <- sd(real_bs_table[,1])
sd_mu_c <- sd(real_bs_table[,2])
sd_sigma_a <- sd(real_bs_table[,3])
sd_sigma_c <- sd(real_bs_table[,4])
sd_rho <- sd(real_bs_table[,5])
sd_psi <- sd(real_bs_table[,6])
sd_vector <- c(sd_mu_a, sd_mu_c, sd_sigma_a, sd_sigma_c, sd_rho, sd_psi)
q_mu_a <- quantile(real_bs_table[,1], c(0.025, 0.975))
q_mu_c <- quantile(real_bs_table[,2], c(0.025, 0.975))
q_sigma_a <- quantile(real_bs_table[,3], c(0.025, 0.975))
q_sigma_c <- quantile(real_bs_table[,4], c(0.025, 0.975))
q_rho <- quantile(real_bs_table[,5], c(0.025, 0.975))
q_psi <- quantile(real_bs_table[,6], c(0.025, 0.975))
q_vector <- rbind(q_mu_a, q_mu_c, q_sigma_a, q_sigma_c, q_rho, q_psi)

mean_mu_a <- mean(real_bs_table[,1])
mean_mu_c <- mean(real_bs_table[,2])
mean_sigma_a <- mean(real_bs_table[,3])
mean_sigma_c <- mean(real_bs_table[,4])
mean_rho <- mean(real_bs_table[,5])
mean_psi <- mean(real_bs_table[,6])
mean_vector <- c(mean_mu_a, mean_mu_c, mean_sigma_a, mean_sigma_c, mean_rho, mean_psi)

med_mu_a <- median(real_bs_table[,1])
med_mu_c <- median(real_bs_table[,2])
med_sigma_a <- median(real_bs_table[,3])
med_sigma_c <- median(real_bs_table[,4])
med_rho <- median(real_bs_table[,5])
med_psi <- median(real_bs_table[,6])
median_vector <- c(med_mu_a, med_mu_c, med_sigma_a, med_sigma_c, med_rho, med_psi)
bootstrap_stats <- cbind(sd_vector, mean_vector, median_vector, q_vector)

write.csv(bootstrap_stats, file = paste0(RESULTS, "/bootstrap_stats.csv"), na = ".", row.names=TRUE)

print("End time")
end_time <- Sys.time()
print(end_time)
time_diff <- end_time - start_time
print(time_diff)
print("finished")

#### EXIT ####

sink()